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Abstract 

We present COSMOCR, a numerical code for the investigation of cosmic ray 
related studies in computational cosmology. The code follows the diffusive shock 
acceleration, the mechanical and radiative energy losses and the spatial transport 
of the supra-thermal particles in cosmic environment. Primary cosmic ray electrons 
and ions are injected at shocks according to the thermal leakage prescription. Sec- 
ondary electrons are continuously injected as a results of p-p inelastic collisions of 
primary cosmic ray ions and thermal background nuclei. The code consists of a 
conservative, finite volume method with a power-law sub-grid model in momentum 
space. Two slightly different schemes are implemented depending on the stiffness of 
the cooling terms. Comparisons of numerical results with analytical solution for a 
number of tests of direct interest show remarkable performance of the present code. 
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1 Introduction 



Structure formation is a major research area in cosmology and a powerful 
observational constraint to discriminate among the viable cosmological models 
[1]. For example, the evolution of the larger structures such as Galaxy Clusters, 
being very sensitive to the underlying mass content of the Universe, is used 
to evaluate the density parameter Q m for matter [2,3]. Numerical simulations 
have proven an invaluable tool and have qualitatively confirmed our theoretical 
framework for the mechanism of structure formation [4]. Today, in the era of 
high precision cosmology, a substantial part of the efforts are directed toward 
the development of a coherent and quantitative picture. This must be inclusive 
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of the feedback from forming structure (e.g., stars, galaxies, and active galactic 
nuclei), which strongly affects the evolution of the observable universe. 

In this context, cosmic-ray (CR hereafter) pressure has been recently recog- 
nized as a possible source of significant dynamical feedback [5] . Indeed, during 
the hierarchical process of structure formation, supersonic gas in- fall and merg- 
ing events invariably generate powerful, large and long-lived shock waves [6-9]. 
These should produce copious amounts of CRs, by way of diffusive shock ac- 
celeration [10], including both electrons and ions. In addition, the post-shock 
gas and diffusively trapped CRs are mostly advected into non-expanding re- 
gions, such as filaments and clusters. It turns out that the energy of most of 
the CR-protons is only marginally affected by radiative losses during a Hubble 
time. It is possible, then, that the latter might accumulate inside large forming 
structures, storing up a substantial fraction of the total pressure there [11]. 
In addition to cosmic shocks other sources of CRs are also possible. These 
include active galactic nuclei, supernovae and stellar winds, all of which are 
candidate for important contributions to the total population of CRs in cosmic 
structures [12]. There is growing observational evidence that significant non- 
thermal activity takes place in clusters of galaxies. This evidence is provided 
by extended sources of radio emission, namely radio halos and relics. From its 
spectral properties and, sometimes, polarization signatures the radio emission 
is interpreted as synchrotron radiation, implying the presence of relativistic 
cosmic-ray electrons and magnetic fields [13-17]. Also, there have been claims 
of detection of radiation in excess to what is expected from the hot, thermal 
X-ray emitting intra cluster medium, both in the hard X-ray band above ~ 10 
KeV [18,19] and perhaps even in the extreme ultra-violet [20-26]. The impor- 
tance of this non-thermal component and its cosmological implications will be 
discussed elsewhere [27,28]. 

In this paper we describe a numerical code that allows a treatment of the 
evolution of various CR populations (i.e., protons as well as primary and 
secondary electrons) in computational cosmology. This code is instrumental 
for a detailed investigation of the aforementioned issues connected with the 
non-thermal activity associated with the formation of the large scale structure. 
In §2 we outline the difficulty of the numerical treatment of CRs in this context 
and a strategy for useful applications. The code is extensively described in §3. 
There we present the algorithm for the treatment of CR ions and electrons 
including their advection in both momentum (§§3.1, 3.2) and physical space 
(§3.4). In addition we describe the inclusion of two source mechanisms, namely 
injection at shocks (§3.5) and production of secondary electrons (§3.6). Finally, 
in §4 we present some numerical experiments to test the code performance. 
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2 Challenges and Approach 



The evolution of the supra-thermal particles, i.e., CR electrons and ions is 
described by the diffusion-convection equation [29]. In comoving coordinates 
the latter takes the form 



at a a z 

bi(p) I + D p ^j- 

+j( x ,P)- 




+ 



where the gradient is with respect to comoving coordinates, x, a(t) is the ex- 
pansion parameter of the universe (such that r = ax is the physical length), 
u = a(t)± is the peculiar velocity (i.e., not inclusive of the Hubble expansion) 
and /(x,p) is the isotropic part of the particle distribution function in co- 
moving units. Here n(p) and D p (p) are the diffusion coefficients in comoving 
coordinates and momentum space, respectively. We recall that on the right- 
hand side of the above equation, the first line contains the spatial terms of 
advection and diffusion respectively; the second line includes the adiabatic 
losses due to cosmic (oc d) and peculiar (oc V • u) expansion, other mechan- 
ical and radiative losses [oc be(p); see eq. (3.7)] and the second-order Fermi 
mechanism (oc D p ); and the third line [j'(x, p)\ represent the comoving source 
term, which accounts for either fresh CRs injection at shock or secondary 
production. Eq. (1) holds only for sub-horizon scales and, more importantly, 
in the non-relativistic limit. These approximation are completely satisfactory 
for the investigations of the problems outlined in the introduction. On the 
other hand, shock acceleration in relativistic flows such as those occurring in 
relativistic jets, gamma-ray bursts and the like, needs to be treated differently 
[30-32]. 

The acceleration and transport of CR electrons and ions, described by the 
above equation, involve physical scales that must be carefully considered when 
a numerical approach is undertaken. The major difficulty arises in the attempt 
to model the behavior of supra-thermal particles nearby shocks. There, the 
smallest length scale of interest is that of the shock thickness, £ s , which, for 
a collision-less shock, is typically of order a few times the Larmor radius of a 
thermal proton [33,34], i.e., 

i a „ TL „J??-~ ( -L-) ( -JL \ 1 io 11 cm. (2) 
ZeB U0 8 K7 \0.1fiGJ V ; 

Here T is the post-shock temperature and B is the magnetic field strength. 
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Particles injected at shocks are thought to be pulled out from the high energy 
tail of the thermal pool and have energy and a Larmor radius a few times 
above the thermal values [35-39]. At even higher energy, CRs have a Larmor 
radius much larger than the shock thickness and, therefore, unlike the ther- 
mal particles, are unaffected by the local shock transition. In the vicinity of 
a shock the dynamics of these supra-thermal particles is dominated by the 
effects of advection and diffusion, which determine the energy gain and es- 
cape probability for the particles [40]. In particular, the spatial distribution 
is characterized by a diffusive scale-length A^(p) = k{j>)/u s (where u s is the 
shock velocity), which determines the distance upstream of the shock to which 
a particle can propagate by diffusing against the advective flow. For a correct 
numerical solution of eq. (1), the hierarchy of the aforementioned scale lengths 
for all relevant momentum, namely that £ s < \d(p)^,^d(Pmax), must be prop- 
erly reproduced in the numerical grid [41]. This demands a full resolution of 
the diffusion length, typically Ax<\d/10, as the shock discontinuity is usu- 
ally spread over 2-3 grid zones in a numerical simulation. When this spatial 
resolution is achieved, the physical roles of diffusion and advection in eq. (1) 
are well accounted for, and the property of continuity of the complete distri- 
bution function, at the shock surface can be reproduced [40]. Notice that the 
continuity of the distribution function across the shock (in the shock frame), 
justified by the fact that the supra-thermal particles do not experience the 
discontinuity, is the very condition that determines the power law behavior of 
the solution /(x,p) emerging from a shock wave [40]. Otherwise, the finite, 
numerical size of the shock, Ax s , induces an error in the solution of eq. (1) of 
order Ax s /\ d [41]. 

However, resolving spatial scales down to the sub-diffusion length of the supra- 
thermal particles turns out to be an extremely demanding task. In fact, if we 
assume for simplicity the case of Bohm diffusion, k = Kb = ^tlv, with tl and 
v ~ c, the particle Larmor radius and velocity, respectively, we find 

where E ~ pc is the particle energy (in the relativistic limit) and u s is the 
shock speed. This means that for electrons with energy ^1 — 10 GeV, which are 
of interest for the non-thermal emission in galaxy clusters, the diffusion length 
is <^ 10 _2 pc and should be even smaller at injection energies. Considered the 
importance of including the cosmic structure on a scale > 50/i _1 Mpc [42], the 
spatial scales to be resolved extend over more than 10 orders of magnitude! 
Therefore, we can typically only afford grid sizes Ax S> A^(p) and following 
the full dynamics of the CRs injection and acceleration at cosmological shocks 
becomes impractical. However, although limited by the above restrictions, one 
can still attempt to find a useful numerical treatment of the CR particles. In 
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particular, the instance that \d{p) / Ax <C 1 also implies that the diffusion 
time at shocks, Td(p) = ^d/ u s < (Xd{p) / Ax) At <C At, i.e., much smaller than 
the computational time-step [41]. The diffusive shock acceleration time, r acc , 
for a strong shock (r = U\ju 2 — 4) is 



r^(p) = — (^ + ^)=3r ( r -±l) # K 20r i(p) 

where r is the shock compression ratio, subscripts 1 and 2 indicate values 
upstream and downstream of the shock, respectively, and for simplicity we 
have assumed K\ = k 2 . Since the time-step for a cosmological simulation is of 
order of 10 8 yr, with the above representative values of magnetic field strength 
and shock speed, At ^> r acc up to a CR energy E ~ 10 15 eV. This energy is 
well beyond the range of energies of interest for CR electrons in clusters of 
galaxies and should include most of the CR ions dynamically relevant in those 
environments. In fact at energies much larger than that, both energy losses 
[43] and diffusive escape [11] become important reducing the population of 
CRs with i?>10 15 eV. Thus, for both CR electrons and ions within the energy 
range derived above, we will assume that shock acceleration is instantaneous 
[36,44] and that a downstream steady-state solution is generated in accord 
with the properties of the shock [41]. The simplest way to implement this 
prescription {acceleration) is to "inject" in correspondence of shock fronts a 
test particle distribution of CRs (cf. [40]) 

f( P ) = finj (-?-) 9S (5) 
\Pinj / 



where the slope is given by 

3r 



(6) 
r — 1 



Here the injection momentum pi n j is determined by the post-shock tem- 
perature (see in §3.5 for details.) In addition, when a pre-existing popula- 
tion of CRs, foip), encounters a shock, the particle distribution is modi- 
fied (re- acceleration) only if the log-slope of the current distribution, q = 
— dlia fo/dlnp, is larger than that determined by the shock compression ratio 
[eq. (6)] (cf. [40]). In general, the injected and re-accelerated particle distribu- 
tions could also be affected by the nonlinear feed-back of the CRs themselves 
[45] . In this case the distribution function deviates from the simple power-law 
given in eq.(5). 
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Once the supra-thermal particles have been produced at shocks, they are car- 
ried through "smooth" flows together with the background gas (before another 
encounter with a shock). The circumstance that A d <C Ax suggests that over 
the scale Ax of the numerical grid, the propagation of the particles is dom- 
inated by the advection term. In fact, the relative contribution of advection 
and diffusion in eq. (1) is of order k/u s Ax = \d/Ax <C 1. Thus, for the spa- 
tial transport of the CR particles in shock-less regions, the diffusion term is 
negligible and can be dismissed in eq. (1) [41]. 



3 The Cosmic Rays Code 

3.1 Advection in Momentum Space 

The difficulty of a numerical treatment of eq. (1) is aggravated by the fact in 
practice we are allowed to use only a "handful" of grid zones in momentum 
space. To illustrate, in the hypothesis that the number of grid zones used for 
the one-dimension momentum variable is the same as for the spatial variables, 
a N A single precision array for a CR population would require already about 
17Gb for N = 256. The number of operations to be computed at each cycle 
grows consequently (oc N 4 ) which is paid at the high cost of a low speed 
code. Finally, we must not forget that the data-set output by the simulation 
becomes huge, rendering the post-computation analysis memory intensive and 
difficult to handle even on well equipped workstations. 

With only a few grid zones spanning the momentum space, a sub-grid model 
of the distribution function is necessary. For this purpose, in the following we 
adopt the general approach first developed in ref. [41], although some aspects 
of the implementation of the numerical schemes below are different from those 
authors. In particular we compute the flux in momentum space exactly up to 
a logarithmic term and provide an alternative scheme for the stiff losses case. 
We sub-divide the region of p-space of interest into N p logarithmically spaced 
momentum bins, bounded by p , . . . ,Pn p - The width of the bins is expressed 
on a log scale as Alogp = Awi = log(pj/pi_i), which is convenient (although 
not necessary) to take as constant. On this grid, we approximate f{p) with 
the following piece-wise power law: 



where and qi are the normalization and logarithmic slope for the % t h mo- 
mentum bin. This approximation is valid when q{p) changes slowly inside each 




(7) 
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bin, which holds true as long as energy losses do not produce strong curvature 
in the distribution function f(p). Then, integration over the i t h momentum 
bin of eq. (1) multiplied by a factor 47rp 2 , gives 



drij 
dt 



where 



Hi 



— V • (urn) + ^ V((/c)Vni) + 

+ l{ { h + 1 v ' u ) p + a ( be ^ + Dp 



dlogf 
dp 



Pi (J^f^-l 

/ 4np 2 f(p) dp = 4tt ft pU ™=2± 

J O C[i 

i-i 

Ipir P 2 ^f d P 

gLp^fdp 



( K )i = 



Pi-i 



Pi 

Qi= J 4vrp 2 j(p) dp; 



Pi-i 



Wf(p) + 

J Pi-1 

(8) 



(9) 
(10) 

(11) 



and where we have used eq. (7) to derive the expression (9). As already pointed 
out, away from shocks we can neglect the diffusive term. Furthermore, we 
henceforth drop the second order Fermi term and focus on the solution of the 
equation 



^l = -V-(un l ) + \b(p)4np 2 f(p)] P ' + Q t 
or L ] Pi-i 



(12) 



where r — t/a and we have introduced 



b(p) = -(^] =(a + Vu) p+alMp) (13) 



dr , . . 

/ tot 



which includes, in addition to other energy loss terms {bf(p); see §3.7), the 
adiabatic terms due to cosmic expansion and fluid divergence (convergence). 
For the sake of clarity, we shall address now separately the details regarding 
the integration of each term of eq. (12). To begin with, we will focus on the 
contribution due to advection in momentum space responsible for the change 



drij 
dr 



b(p) 4np 2 f(p) 



Pi 



Pi-i 



(14) 
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Time integration of eq. (14) yields 

< + At -< = -At (d>? - (15) 



where 

r+Ar 

= | fe(p)47rp 2 /(r , ,p)Urfr', (16) 

T 



is the time- averaged flux and the integrand is evaluated at the cell boundary 
Pi. Recalling eq. (7), the above integration is readily carried out after using 
eq. (13) to rewrite it as 

Pu 

^ = £1 p 2 fM dp (17) 

Pi 



where 



J = < 



i + 1 if b(pi) > 
% if b{pi) < 



(18) 



and p u is the upstream momentum, solution of the integral equation 

Pu 



After updating n« based on (15), and including the contributions from spatial 
advection and injections (detailed below), the sub-grid model of the distribu- 
tion function (7) is reconstructed by solving for fi and qi at each computational 
cell. For each new value n^, fi and qi are related by eq. (9), so one additional 
constraints is necessary. For the second relation we assume that the curvature 
of the spectrum is constant, i.e., q i+ i — qi = ^ — (cf. [41] for more detail on 
this). This is a sensible approximation when cooling is weak. It is, therefore, 
very suitable for cosmic ray ions, which, in a cosmological context and for the 
energy range of interest here, suffer only minor losses even during a Hubble 
time. However, for the case of electrons, which cool much more rapidly, this 
approximation breaks down and additional measures must be taken, which we 
address below. 
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3.2 Alternative Scheme for Fast Cooling Electrons 



The assumption of a smooth CR energy distribution (constant spectrum cur- 
vature), which allows one to infer the slope of the distribution function from 
the particle number density in each bin, breaks down in the presence of strong 
cooling. This is an important issue for electrons, since the cooling times as- 
sociated with inverse Compton and synchrotron losses are much shorter than 
the time-scale of a cosmological simulation. It is important to notice that the 
effect of strong inverse Compton or synchrotron cooling is that of producing a 
cut-off in the particle distribution, not just a mere steepening of it. Since each 
momentum bin spans a considerable range of values, it is often the case that 
the cut-off falls in the middle of a bin, so that only the part of the bin below 
the energy cut-off is populated, while upper part is completely depleted. This 
situation is likely to occur for any of the momentum bins, not just those at 
the highest energies, because, due to additional Coulomb losses (see §3.7), the 
cooling times for the electrons at all energies of interest are always smaller 
than or comparable to the Hubble time. For the same reason, except around 
shocks where they are continually injected, primary CR electrons in most of 
the simulated volume only occupy a very few of the lower energy momentum 
bins. This adds a non-trivial complication because at least three momentum 
bins are necessary for the reconstruction scheme mentioned in the previous 
section to work [41]. 

In order to circumvent these difficulties, we have devised an alternative scheme 
to be employed for the electrons only. The idea consists of replacing the "con- 
stant curvature" assumption with a different condition that allows us to recon- 
struct both normalization, f\, and slope, qi, of the piecewise power-law distri- 
bution function (7) in each bin. A natural constraint is provided by the physical 
condition that the total energy be conserved, except for the explicit sources 
and sinks. The corresponding equation is derived by taking one moment of the 
transport equation (1). Thus, after multiplying by a factor 4irp 2 T(p), where 

T{p) = (7 — \)m e c 2 is the particle kinetic energy (where 7 = l/^/l — (v/c) 2 
is Lorentz factor), and integrating over the i t h momentum bin, we have 



dej 
dt 



-- V-(u £i ) + 4 V((K T )iVei) 

1 

a 



a+ - V • u ) p +abi(p) 



Wf(p)T(p) 
4ttp 2 /(p) 



Pi-i 
p 



Pi-i 



\jm\c 2 + p 5 



dp 



+Si. 



(20) 
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As before, we have dropped the second order Fermi term, and introduced the 
following definitions 



J 4 Qi 



Pi-i 



(21) 



{ Th ~ WyNWrndTp (22) 

Si= J 4vrp 2 j(p)T(p)dp, (23) 

Pi-i 

where eq. (7) has been used and the sub-relativistic contribution has been 
ignored in order to derive the last expression in eq. (21). The first two terms 
on the right hand side of eq. (20) represent the usual spatial advection and 
diffusion respectively. The third and fifth terms account for advection in mo- 
mentum space and external energy source respectively, analogous to the cor- 
responding third and fourth terms in eq. (8). Finally, in the fourth term on 
the right hand side we combine, for reasons of numerical convenience, both 
the CR pressure and the sink contributions. 

As before, the diffusion term can be neglected, whereas both spatial advection 
and energy injection will be treated in the next sections. The changes affecting 
the distribution in momentum space are then provided by the third and fourth 
terms according to 

d£ i \u \ a 2 ,/ \ rrii slK ? u \ 4vrp 3 /(p) 



Or 



b(p)47rp 2 f(p)T(p)] Pi f b(p) dp, (24) 

" '' ■' ■ ID; c' 2 ■ ji 1 



Pi-i V e 



where r = t/a and b(p) is defined by (13). For the range of energy of interest 
here p ^> m e c. Then, by using the prescription in (7) for f(p), the second 
term in eq. (24) can be rewritten as EiRi(qi,pi-i), where: 

Pi 3 _ q _ Pi 

Ri(qi,Pi-i)= I b{p) j2== dp I f p 2 -i>T(p)dp. (25) 
P i-i ^/m 2 e c 2 +p 2 v l x 



Thus, after integration over a time-step eq. (24) reads 

eJ +AT (l + - eT (l - ^ ) - Ar ($f - (26) 
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where the term oc Ri has been integrated implicitly and 



i- J b(p)4np 2 f t (r',p)T(p) \ Pi dr', 



t+At 




(27) 



r 



is the time-averaged flux whose integrand part is evaluated at the cell bound- 
ary pi. In analogy with the previous section, by using the definition (7) and 
eq. (13) the above integral is readily rewritten as 



where the index j and the upper integration limit, p u , are defined in (18) 
and (19) respectively. In addition, the maximum momentum of the electrons 
distribution, p cut , is followed explicitly by means of equation (13), which can 
be easily integrated analytically [46] for the cooling mechanisms of interest 
here. Therefore, in this scheme in eq. (28) and eq. (17) we adopt the minimum 
between p u and p cut . Following explicitly the value of p cut , i.e., the high energy 
cut-off of the distribution, allows us to describe the electrons in the last mo- 
mentum bin by a normalization and slope appropriate for the still populated 
portion of that bin. 

After accounting for the advection in physical space and source terms treated 
in the following sections, we have information on both the particle number 
density and kinetic energy at each momentum bin. Taking the ratio of energy 
to number density for each bin, we find for pi ^> mc. 



which can be solved in the unknown qi through Newton-like method with only 
a few iterations. The value of fi, which is actually not explicitly needed in any 
integration step, can easily be computed from either definition (9) or (21). In 
this method both fi and qi are derived from quantities pertaining exclusively 
to the i t h bin and no information about the adjacent bins is required. Thus 
the scheme works with even one single cell or a fraction of it (since we keep 
track of p cut which can even be smaller than pi). This allows us to follow for 
cosmological times the low energy electrons with 7 ~ 300. The electrons in this 
energy range have the longest lifetime against energy losses in a cosmological 
environment [47]. Thus, about 10 9 yr after a CR distribution has been produced 
in a shock, these are the only electrons left. However, It is important to keep 
track of them, even if most of the electrons at higher or lower energy have been 






(29) 
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depleted, because they can contribute to the extreme ultra-violet radiation 
detected in clusters of galaxies. 

3.3 Constraint on the Time-step 

In order for the scheme described in §3.1 and §3.2 to behave properly the 
time step At must be such that I log — I < ecFL,Aw where ecFL < 1 is a 

Pi 

Courant-like number. Given the very long cooling times for the proton CRs 
compared to the typical cosmological time-step, there is no need to enforce 
the above condition for this case. For the electrons, however, we find that in 
order to follow accurately the evolution near the high energy cut-off a limit 
on the time-step At <; 0.lT sync+IC should be adopted. 

3.4 Advection in Physical Space 

While we have focused so far almost exclusively on the integration of the 
diffusion-convection equation in momentum space, the full evolution of the CR 
particles must also include the effects of spatial transport. Since we neglect the 
spatial diffusion for reasons already pointed out in §2, the main contribution 
in this respect is provided by the advection terms. In eq. (12) and (20), these 
are of the form 



with Q = rij or £j respectively. These terms are integrated by means of a van 
Leer method [48], which is a conservative, second order accurate and Total 
Variation Diminishing scheme. The scheme allows integration of the advec- 
tion terms in one dimension and the full integration along all three different 
directions is achieved by directional splitting [49] . As a general feature of con- 
servative, Godunov-like schemes [50], the quantity that is being updated is 
the volume-average over the computational cell Q , i.e., 



The time update is accomplished through a double integration of eq. (30) over 
the cell volume and over the time-step At. Upon performing such integrations, 

1 We point out, for the sake of completeness, that because of this volume- averaging 
step, the distribution function / employed in §3.2 should be read as /, i.e., its 
volume average. 



dcj 
8t 



V-(uCi) 



(30) 




(31) 
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we find 



-At £ (32) 

.£=1,2,3 



where the index I runs over the three spatial coordinates and accounts for 
the flux through the computational cell boundaries along the three spatial 
directions. The time-averaged flux in physical space is defined as 



t+At 



t 



with the integrand evaluated at the cell boundaries X£. The numerical solu- 
tion c(t, x,p) is reconstructed from the volume average values (c) as a piece- 
wise linear interpolation. An important part of this step is the effort to re- 
strict the slopes of the linear interpolation. In fact, having the potential to 
become artificially large near extrema or discontinuities, they can cause un- 
wanted oscillations and render the scheme unstable. This problem is elim- 
inated by demanding that the interpolated function is monotone [48]. The 
reconstructed solution can then be upgraded "exactly" in Lagrangian coordi- 
nates as c(t',x,p) = c(t, x — \l vdt,p), allowing the calculation of the flux in 
eq. (33) at each interface xe. 

3. 5 Injection 

The source terms in eq. (12) and (20) are responsible for the variations 

n^-nl = Y t (34) 

gr+Ar _ E r = (35) 

respectively, where 



^ t+At t+At 

Y '=vKt I I Qid3xdt= Ai I Q ' dt (36) 

t V t 
t+At t+At 

' ' r Sid 3 xdt = — [ Stdt (37) 



VAt J J At 

t v t 



(38) 



and Qi and Si have been defined by eq. (11) and (23), respectively. In this 
section, we shall focus our attention on source contributions provided by CR 



13 



injection at shocks. With this term we indicate the number of particles that 
upon passing through a shock are assumed to undergo the diffusive shock 
acceleration mechanism. Since the latter is treated as instantaneous, injection 
here basically refers to the deposition in the post-shock region of CRs with a 
power-law distribution in momentum space which accords with the diffusive 
shock acceleration theory (see §2) for further details). Next (§3.6), we will 
outline the scheme for the injection term due secondary electrons produced in 
inelastic p-p collisions. 

The scheme described here for the injection of CRs at shocks is based on the 
thermal leakage model which seems to be observationally supported at least 
by in situ measurements at the earth bow shock [35]. Here, we assume that 
upon passing through the shock most of the gas thermalizes to a Maxwellian 
distribution with post-shock temperature T 2 , i.e., 



with n,2 the total number of particles of the distribution. For a particle to re- 
turn upstream it is necessary not only that it propagates faster than the shock 
wave, but also that it has enough energy to escape "trapping" by Alfven waves 
generated in the downstream turbulence [51,52]. Thus, only those particles in 
the high energy tail of the thermal distribution will have a chance to re-cross 
the shock and get injected into the acceleration mechanism. The numerous, 
complicated details of the physics underlying the injection mechanism are con- 
veniently modeled by a few parameters [37-39,53]. One of them, ci, defines 
the momentum threshold for the particles of the thermal distribution to be 
injected, as 



In practice we assume that at p in j the thermal distribution, f m (pinj) [ ec i- (39)], 
joins "smoothly" into the power-law distribution of the CRs. That implies 



where q = 3r/(r — 1), r is the shock compression ratio and pi n j < p < p max . 
With these choices, C\ is the only free parameter in the injection model. In 
fact, the total number of injected particles is given by 




(39) 




(40) 




(41) 




(42) 



Pinj 
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and the injection efficiency, i.e., the fraction of downstream thermal gas par- 
ticles that are injected in the acceleration process, is [cf. eq. (39) and (40)] 



Vinj 



n 2 



I Pmax 
\ Pinj 



max ( 1 



(43) 



practically independent of T 2 . In an alternative approach, where the down- 
stream gas distribution is not specified, Vinj and C\ can be taken as indepen- 
dent parameters [38,45]. With the above setting, CR particles are injected per 
unit shock surface at a rate 

- _ pi u s _ n 2 u 2 , s 

Vi Vinj Vinj \ ) 

m p r 



where n 2 u 2 = ri\ U\ is the number flux of thermal particles impinging on 
the shock. Given the injection rate, Q iy and the distribution of the injected 
particles in eq. (41), the energy injection term, Si, can be easily computed as 
well by means of eq. (23). 

The actual values of c\ and r] in j are not known and most likely they are 
not constant, but a function of the conditions of the flow. For ionic CRs, 
observations at the earth's bow shock indicate that an injection efficiency, rjinj, 
can be as high as 1.25 x 10~ 3 [54] or even ~ 10~ 2 [35]. In theoretical studies of 
shock acceleration at supernova remnants the parameter C\ is assumed in the 
range 2.3-2.5 corresponding to values of rjinj ranging between a few xlCT 3 to 
10~ 4 [55,37,53,38,51,39,45]. It is important to notice that the value of c\ (in 
addition to p max and the slope q) , regulates the amount of flow kinetic energy 
that is transferred to the CRs. CR protons can be dynamically important by 
exerting a pressure 

Pmax Pmax < 

4"7T f Atc f V 

PcT = T I m p3 '" ip = T C / m (m^+p^ *' (45) 

Pinj Pinj 



The ratio of the CR pressure to the ram pressure of the flow, p\u\, can be 
regarded as a first order indication of the relative importance of the two com- 
ponents. For a flat CR distribution, i.e., q ~ 4, which is typical of the cosmo- 
logical case [5,28], by using the above expression for P cr and neglecting the 
sub-relativistic contribution we find 



\ 3— q o I P max 

P cr _ 8 2 3 2c? / m p c\ ( c y { Pm] 




P\U\ 6 V 7T \p inj I \U t 



(46) 



Thus, for a CR ion component extending between 1 GeV and 10 6 GeV, ap- 
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propriate for shocks in galaxy clusters, the back-reaction of the accelerated 
particles would not be negligible with the aforementioned choices of the pa- 
rameter c\. It appears that in cosmo logical simulations of large scale structure 
formation, most of the kinetic energy is processed by shocks with Mach num- 
ber in the range 3-6 [9]. According to recent results from a newly developed 
Adaptive Mesh Refinement scheme with a resolution down to the diffusion 
length scale and capable of including self-consistently the dynamical role of 
the accelerated particles [56] shocks with Mach number M ~ 3— 6 are expected 
to be "moderately" efficient (up to 30%) but not strongly modified. Thus, the 
effect of the back-reaction can be emulated by slightly increasing the value of 
Ci, which has the effect of reducing the total number of injected CRs. For the 
above range of flow parameters, we find that a value of c\ = 2.6 produces an 
injection efficiency r) in j consistent with above nonlinear calculations [56,57]. In 
general, however, one must be very cautious about the choice of C\. Its value 
should be determined on a individual basis, based on the properties of the sim- 
ulated shocks where most of the CRs are being produced and in accord with 
the possible non-linear behavior of the shock acceleration mechanism there. 

The physical process of injection of CR electrons is more complicated and 
basically not yet fully understood. Part of the reason is due to the fact that 
electrons carry only a small fraction of the momentum of the flow and, there- 
fore, they do not affect the dynamics of the shock. Thus, unlike the ions, 
whose behavior can be constrained by general considerations of energy and 
momentum conservation, the electrons behave as test particles and their dy- 
namics is determined by the details of the plasma wave-particle interactions. 
The latter is very difficult to model appropriately in computer simulations. 
In addition, until recently the observational results available in this respect 
were very limited. In this situation, the progress in the understanding of CR 
electron injection at shocks has been relatively slow. However, for the present 
purpose of a numerical treatment of CR injection at cosmic shocks, a viable 
and reasonable approach is to assume that the ratio of the cosmic ray electrons 
to ions at relativistic energies is fixed to a value R e / P [58] . This phenomenolog- 
ical approach is supported by some observational evidence from experiments 
with Galactic cosmic rays indicating that this ratio is possibly in the range 
1-5 % [59,60]. 



3. 6 Production of Secondaries 

Electrons and ions injected at shocks provide the main source of CRs and are 
usually referred to as primary. As high energy protons propagate through the 
galactic or intergalactic medium, they collide with the background thermal 
protons and, from the hadronic interaction, secondary products are gener- 
ated. The latter include photons, leptons and hadrons and, therefore, might 
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be an important source of relativistic electrons. The main channels for the 
production of secondary electrons (positrons) are [61]: 



p + p^n^ + X (47) 

p + p^K^ + X. (48) 

where X indicates all the by products of the reactions. Charged pions decay 
with a lifetime of 2.6 x 1CT 8 s [62], mostly into 



7T + ^/i + + Z/ M (49) 

n'^fi' + v^ (50) 

Kaons, analogously, have a lifetime if 1.24 x 1CT 8 s and decay primarily into 
muons (63.5%) and pions (21.2%) [62]: 



K + ^pL + + ^ (51) 
K'^iT + Vv (52) 
A' ± ->7r + 7r ± (53) 

The relative contributions of the two channels (p + p — > + X and p + p — > 
K ± + X) are energy dependent. So, for example, the fraction of secondary 
muons from K decay is 8% at about 100 GeV, 19% at 1 TeV and approaches 
asymptotically 27% at higher energy [61]. In turn, muons have a lifetime of 
2.2 x 10~ 6 s [62] before decaying into 



+ 



+ v e + Up 



e +v e + v, 



(54) 
(55) 



In addition to p + p inelastic collisions, the above cascades are also triggered 
by the interaction of p+He, a+H and a+He, which, for example, increase the 
overall yield of secondary e ± by a factor 1.4 for a metal composition relative 
to the interstellar medium [63] . 

In general we can write the production spectrum of secondary electrons as [64] 

-max/- \ 

j s (e s )=n H J de p J p (e p ) (Ccri(e p )) J dSi F s (e s , e i )F i (e i , e p )(56) 

i 7T)K ffiiji -rain I - \ 

b p £ i \ £ s) 
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where J p (s p ) is the proton flux; ((<Ji(e p )) is the inclusive cross section^ of 
the processes (47) and (48); e p nm is the minimum proton energy required to 
produce a meson of energy e™ m and ef bax the maximal energy of the produced 
meson; e™ n in turn is the minimum required energy of a meson for produc- 
tion of secondaries of energy e s ; finally F^e^k, e p ) are the spectra of n and 
K produced from the collision of a proton of energy e p and F s (e s , £n,K) the 
distribution of secondaries from the subsequent decay of the above collision 
products. For a full summary of the technique employed here in calculating 
the cross sections and the functions F s and Fi we refer to ref. [64]. 



3. 7 Energy Losses 



We have considered the energy losses due to a variety of physical mecha- 
nisms that become relevant at different energy regimes. For the electrons, the 
most effective process is due to Coulomb losses in the low energy end and 
synchrotron and inverse Compton emission at high energies. Bremsstrahlung 
losses are also included for completeness, although less relevant. For ion CRs, 
Coulomb losses, which are mechanical in nature, are dominant below relativis- 
tic energy. However, given the much smaller cross section for radiative losses 
for the ions (oc l/m^) as compared to electrons (oc l/m 2 ), the next important 
loss mechanism for the CR ions beside adiabatic expansion is that due to in- 
elastic collisions with the thermal background nuclei. In the following we shall 
provide the functions 

b(p) = f t (58) 



and the associated cooling time 



T~cool — 7 / \ (59) 

b(p) 



relative to each relevant process, and separately for electrons and ions. 
3. 7. 1 Electrons 

Losses due to Coulomb collisions are expressed by [65] 



2 Inclusive means that which describes the process 

p + p^i + X (57) 
where i is in general a secondary particle 
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3 



\ dt Jcoul m -° 2 

= 3.01 x 10~ 29 jl + [ln(l + p 2 ) 1/2 - lnn ]^5g} nergcm-(60) 

where Z is the electric charge of the ion species, n is the number density of 
the background gas in cm -3 and p = p/m e c. Taking Z — 1, the corresponding 
cooling time is defined as 



Tcoui = 2.53 x 10 £ 



\W0J Vl0~ 3 cm- 3 



yr. 



(61) 



Bremsstrahlung losses are defined as [66,65] 



dp^ 
.dt. 



Aa.fr 2 'yrrieC 2 Z[Z + 1) n 



brem 



ln(2 7 ) - - 



3.8 x 10~ 33 {ln[2(l+j3 2 ) 1/2 ] - i| pnergcm" 1 



(62) 



where Oif is the fine structure constant, r e the classical electron radius. Such 
mechanism is basically unimportant as it becomes effective on a time-scale 



Ur 



6 x 10 



12 



n 



10" 3 



cm" 



(63) 



again taking Z — 1. The other relevant loss mechanisms for electrons are 
due to synchrotron emission and inverse Compton scattering off the cosmic 
microwave background photons. After averaging over the electrons pitch angle, 
their combined contribution is given by 



= 8.94 x 10' 25 {u B + u cmb )p 2 erg cm" 1 (64) 

where U B and u cmb = 4.2 x 10~ 13 (1 + z) A are the energy density in magnetic 
field and cosmic microwave background at a given cosmological red-shift, z, 
respectively, both in units of erg cm" 3 . For high energy electrons this is the 
most severe energy loss mechanism with a typical time-scale at red-shift z = 

« = 2.3X10' (±) '(l + ^f-)" yr- (65) 
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3.7.2 Ions 

Coulomb collisions are efficient at a rate [65] 



'dp\ 2nZ 2 e i , / 7 2 m=!c 4 

— = — n In 1 

dt) 
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p 4 /(l+p 2 ) 



1 + 



ln 



2(m e /m p )(l+p 2 y/ 2 / 



— Inn 



n erg cm 



-i 



75.7 



(66) 



where beta = v/c, 



X r 



3tt\ 1/3 / 2kT e 
4 / \ m P c 2 



1/2 



64 x 10 " ( Jk) 



(67) 



and we now define p = p/m p c. For protons in the low energy end within the 
range considered here, Coulomb collisions are the only significant energy loss 
mechanism, with a time-scale comparable to the Hubble time, 




n ' 1 



r Coul = 4.8 x 10 9 ( JL ) (^5=3^=3 ) yr, (68) 



where Z — 1 is used. Photo-pair and photo-hadron production of CR ions A, 
interacting with the cosmic microwave background photons, i.e., 

A + 7 cm6 ^A + e + + e~ (69) 
A + 7cm 6 -> A + yr (70) 

are mostly negligible at the energies we include. In fact, for red-shift z = and 
taking the average CMB photon energy (e cm b) = 7 x 10~ 4 eV, the thresholds 
for the above reactions, in terms of the Lorentz 7 factor, are respectively 

2 

7mm = T~~~~~T ^ 7 X 10* (71) 

7m ,„ = !p4((l + !^Wxl0» (72) 

reachable only for ultra high energy CRs (where £ is the multiplicity of the 
produced pions). Rather, inelastic collisions of CR ions off nuclei of the thermal 
background gas are more significant, being at a rate [43] 
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= <J n ,inel{E p ~ VTl p (?)n 

p-p 

2.91 x 1(T 29 (p - P)n erg cm" 1 if e v > 1.22GeV 

V 7 p " (73) 

if e„ < 1.22GeV 

V 1 

and characterized by a energy loss time-scale 

T "= 5 - 5xioi °(io^r yr < 74 > 



4 Test Results of the Code 

In the following we present a set of numerical experiments that test the perfor- 
mance of the code in various cases of direct interest. In particular we consider 
the time evolution of an initial distribution of CR protons and electrons as 
they lose energy subject to mechanisms described in §3.7. Alternatively, we 
also present the evolution toward the steady state of electrons which are con- 
tinually injected as secondary products of p-p collisions and at the same time 
lose energy due to Coulomb and inverse Compton losses. For the first experi- 
ment we take the initial spectrum of the cosmic rays (protons and electrons) 
as a power-law with logarithmic slope q = 4.3. A similar distribution (but 
with a different normalization) will also provide the parent CR protons that 
produce the secondary electrons in the last example. The numerical solution, 
which gives the number density of CRs in each momentum bin as a function 
of time, t, is plotted and compared with the analytical solution. The latter is 
obtained by integrating between each momentum bin bounds the exact, ana- 
lytical CR distribution corresponding to the time t [67], for the same initial 
distribution evolved by the code. As for the Coulomb, bremsstrahlung and 
p-p collision losses, we assume that the CR protons or electrons propagate 
through a medium with number density n gas = 10~ 3 cm~ 3 , typical of the core 
of clusters of galaxies. In order to further mimic cosmic conditions, the test 
are run for a time t h = 1.5 x 10 10 yr so that all the cooling time-scales of 
relevance in a cosmological simulation are included. 

4-1 Evolution of Cosmic Ray Power-Law Distributions 

In fig. 1 we present the temporal evolution of a distribution of CR pro- 
tons subject to Coulomb losses and p-p inelastic collisions. For this experi- 
ment we use 8 momentum bins, with the CR spectrum extending between 
Po = O.lGeV/c up to p 8 = 10 6 GeV/c (c is the speed of light). The dis- 
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Fig. 1. Evolution of CR proton spectrum with initial logarithmic slope go = 4.3, 
subject to Coulomb losses and p-p inelastic collisions over a period of 15 Gyr. 



tance between subsequent bin bounds is constant on a logarithmic scale and is 
Aw = [log(pg) — log(po)]/8 = 0.875. Since the losses are only mild in this case, 
the numerical evolution of the CRs is computed by the scheme presented in 
§3.1. For each momentum bin, on the ordinate axis we plot the number density 
of CR protons as computed from the code (pentagons) and from the analyt- 
ical solution (square). Here and in the following figures, in order to facilitate 
the visual comparison of the results, for each bin the squares are drawn in 
correspondence of the middle of the bin and pentagons at a slightly higher 
momentum value, although the bins are identical for both cases. The points 
in the plot correspond to the following times t/r H = 0.1,0.3,0.5,0.8,1. The 
time sequence is such that, for each bin, points in the upper positions corre- 
spond to earlier times. As we can see, above 1 GeV the CRs are only sensitive 
to the p-p losses which deplete all the momentum bins at the same rate. In 
this regime the two solutions, the numerical and the analytical one, are, for 
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Fig. 2. Evolution of CR electron spectrum subject to Coulomb, bremsstrahlung and 
inverse Compton losses over a period of 15 Gyr. The initial logarithmic slope was 
set to 4.3 

all practical purposes, identical. At lower energies, i.e., for the first bin, the 
numerical solution still follows very accurately the exact solution. Only in the 
last two points of the time sequence, after t — 1.2 x 10 10 yr = 0.8 t Cou i (p/0.3), 
the numerical solution cools a little bit faster than the analytical one and the 
number of CR protons in the first bin (only!) is slightly underestimated by a 
factor ^ 1.4. 

Analogously, in Fig. 2, we present a comparison between the numerical (pen- 
tagons) and the analytical (square) solutions for the evolution of an initial 
power-law distribution of CR electrons. In this case, the numerical solution 
of the evolution of the CR distribution is computed by means of the more 
sophisticated (and more expensive) scheme outlined in §3.2. Again we use 
8 momentum bins, but the extrema of the CR electrons spectrum are now 
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Fig. 3. Evolution of secondary CRs produced in p-p collision by a power-law distri- 
bution of parent CR protons with logarithmic slope was set to 4.3 



Po = 50MeV/c and p$ = 2 x 10 5 MeV/c with a distance between subsequent 
bin bounds Aw = 0.45. The inverse Compton losses are computed by setting 
u>b <^ u cmb = 4.2 x 1CT 13 . The values of the CR number density for each bin 
in Fig. 2 has been plotted for the following times (ages) of the distributions: 
t/ TH = 5 x 1(T 4 , 1 x KT 3 , 3 x KT 3 , 1 x 1(T 2 , 2 x KT 2 , 0.1, 0.25, 0.45, 0.6, .7. The 
above times have been chosen in order to test the accuracy of the numerical 
routine in different regimes. For example, the first three times correspond to 
a fraction of r IC + S ync for particles with p ~ 10 4 , whereas the last times are a 
significant fraction of t Cou i- Again, since the evolution of the CR spectrum is 
determined by energy losses, the time evolution is such that the upper points 
correspond to earlier stages. Since the cooling time for the CR electrons in the 
high energy bins is much shorter than the total evolution time th, eventually 
such bins will be evacuated. Since the plotted dynamical range is limited, we 



24 



adopt "upper limit" symbols for those bins where cooling has reduced the 
number of particles below an arbitrary numeric level of 10~ 4 . We notice once 
again a very accurate correspondence between the code and the analytical so- 
lution. As for the CR protons, the numerical solution departs slightly from the 
analytical one at the lower energies (for the first 2 bins). However, this now 
takes place only after t > 0.6 th — 7.1 rcouh to be compared with t > 0.8tcw 
for the protons case. This means that the present scheme is more accurate not 
only for the evolution of the high energy part of the spectrum but also the 
low energy end. This feature is relevant in the context of cosmological simula- 
tions because the time-scales for the energy losses of electrons at all energies 
of interest are significantly shorter than a Hubble time. Notice that, unlike 
the previous case, the numerical solution now overestimates the total number 
of CRs in the lower momentum bins. The discrepancy, however, appears only 
when a momentum bin has already been significantly evacuated of electrons 
and is close to complete depletion. From Fig. 2 we can see that the difference 
between the true and the numerical solution (for both first and second bins) 
is only a small fraction of the difference between the initial and the current 
value. Also, notice that after the second momentum bin has been evacuated, 
i.e., around t/rn ~ 0.5, the scheme is working with only one momentum bin. 

Finally, in Fig. 3 we compare the numerical and analytical results for the 
evolution of a distribution of secondary CR electrons, continuously injected 
through the processes described in §3.6. In this case we start from a low 
background distribution of CR electrons and the continuous injection of new 
particles builds up the steady state distribution. Thus, unlike in the previous 
examples, now the time evolution is such that the lower points correspond to 
earlier times. The CR electron distributions have been plotted for both the nu- 
merical (pentagons) and analytical (square) solutions for the following times: 
t/ TH = 5 x 10" 4 , 1 x 10" 3 , 3 x 10" 3 , 1 x 10" 2 , 2 x 10" 2 , 0.1, 0.3, 0.5, 0.8, 1. Again 
the numerical scheme follows very precisely the evolution of the spectrum of 
the secondary electrons. Since the slope of the parent CR protons was q = 4.3, 
the injection spectrum has a logarithmic slope qi = 4.3 [43]. Thus, as expected, 
the slope of the steady state spectrum is respectively steeper at high energies 
and flatter at the low energy end than the injection spectrum by one unity. 
We notice, however, that the first and last bins are even flatter and steeper 
respectively than the above prediction. At low energies, this is a physical effect 
due to the reduction of the injection rates of secondary electrons as we ap- 
proach the lower energy limit for electron production ~ ~ 50MeV [68]. 
The steepening of the slope of the last bin, on the other hand, is a boundary 
effect, namely the lack of incoming flux of electrons from higher energies. This 
is due to the fact that we included an injection spectrum only up to electron 
energy of 100 GeV. Since the same injection spectrum was used for both the 
numerical and the analytical solutions, the test comparison is unaffected by 
this choice. 
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Fig. 4. Gray-scale image of a 2-D slice of the gas density distribution in units cm 3 
for a cluster selected from the simulation. Linear size of the image is ca. 5 ft, Mpc. 

4-2 Cosmic Ray Distributions in a Galaxy Cluster from a Cosmological Sim- 
ulation 

At last we present the distribution of CR protons as well as primary and sec- 
ondary electrons, produced in a cosmological simulation of large scale structure 
formation performed by implementing COSMOCR into a hydro+N body cos- 
mological code [69]. The simulations are described in detail elsewhere [5,27,28]. 
For the present purpose it suffices to mention that we assumed a SCDM cos- 
mology with an initial power spectrum of density fluctuations with cluster 
normalization <r 8 = 0.6 and spectral index n — 1; normalized Hubble constant 
h = H /(100 km s _1 Mpc -1 ) = 0.5; total mass density tt M = 1; and baryonic 
fraction Qj, = 0.13. For the simulation we selected a cosmological volume de- 
fined by a cubic region of comoving size 50 ft, _1 Mpc and use 256 3 cells for the 
baryonic matter and 128 3 dark matter particles. This corresponds to a spa- 
tial resolution of ~200/i _1 kpc. The CR protons and primary electrons have 
been injected and accelerated at shocks formed during the simulation with the 
prescription described in §3, whereas secondary electrons are produced in p-p 
collisions. 

We have selected one of the clusters formed in the simulations. In Fig. 4 we 
show a gray-scale image of a slice of the gas density in units cm~ 3 through the 
core of such clusters. For the same cluster, in Fig. 5 (left panel) we present 
the analogous distribution of CR protons, and in Fig. 6 (left panels) that 
of primary (top) and secondary electrons (bottom) as well. In addition to the 
CR spatial distribution for each species we also project onto the corresponding 
right panel the CR spectral distribution (upper-half) and log-slope (lower-half) 
for two locations in the intra-cluster medium: the cluster center (open squares) 
and its periphery nearby an accretion shock (filled pentagons). 

The spatial distribution of each species resembles that of the gas density. In 
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Fig. 5. Left: distribution of the total number density of CR protons in units of 
cm~ 3 . The length of a side of the image is 5 /i _1 Mpc. Right: number density (top) 
and power law index (bottom) for each momentum bin for two locations in the 
intra-cluster medium, at the cluster core (open squares) and periphery (filled pen- 
tagons) . 



Fig. 6. Same as in Fig. 5 but for primary (top panels) and secondary electrons 
(bottom panels). 
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addition, CRs are only present within the shocked regions. This makes sense 
because roughly speaking according to our injection mechanism the CR parti- 
cles are generated downstream of shocks and are then passively advected with 
the flow. For the case of the protons (Fig. 5) the number density is higher in 
the cluster core (right panel - open squares) than in the periphery (right panel 
- filled pentagons) due to the adiabatic compression undergone by the CRs 
together with the gas during the formation of the cluster. And in fact, in the 
inner regions of the cluster where the gas is denser and, therefore, the Coulomb 
losses stronger, the CR proton spectrum shows a mild flattening at low ener- 
gies (Fig. 5, right-panel bottom-section), with respect to the distribution in 
the outer region. 

Analogously, the total number of CR primary electrons is larger at the cluster 
core than at its periphery (Fig. 6, left-top panel). However, at the high energy 
end the primary electrons in the cluster center have been completely depleted 
due to severe inverse Compton losses, unlike in the peripheral region behind 
the accretion shock where fresh CR electrons are being generated. Such differ- 
ence is reflected also in the values of the log-slope of the distributions, steep 
for p > 10 3 in the cluster center and still relatively flat at the periphery near 
the shock (Fig. 6, left-top panel). Finally, the bottom part of Fig. 6 shows the 
properties of CR secondary electrons. Again the number density is higher at 
the center of the cluster than at its outskirts. Notice that the number density 
of secondary electrons roughly scales In fact, this qualitative 

difference is reflected, approximately in the right proportion, in the larger dis- 
tance between the two density curves in the bottom right panel (upper-half) 
of Fig. 6, as compared to those in Fig. 5. Notice also that the log-slope of 
the distribution of the secondary spectrum is steeper (flatter) than that of 
the parent CR protons at high (low) energies by approximately one unit as 
expected. Again, the first and last bin make an exception to this expectation, 
in accord to our explanation in §4.1. 



5 Conclusions 

In this paper we described COSMOCR, a code for cosmic ray associated stud- 
ies in cosmological simulations. We have described the scientific motivations 
behind such code, the challenges posed by the realization of a suitable algo- 
rithm and the various numerical schemes adopted here in order to follow the 
evolution of different CR species. We have shown in §4 that the present code is 
able to follow with high accuracy the evolution of a CR distribution of protons 
as well as primary and secondary electrons subject to the mechanisms of en- 
ergy loss of interest in the intra-cluster medium. Additionally, COSMOCR has 
been implemented in a cosmological code [69]. The resulting distribution of 
CR populations are in accord with what we know about the properties of the 
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cosmic shocks responsible for the acceleration of such CRs [9], and with our 
knowledge of the transport and losses mechanism that regulate their evolution. 
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